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INTRODUCTION 


Reference  1  presents  the  general  theory  of  canonical  eigenvalues.  This  report  applies  that  theory 
to  underwater  sound  propagation  in  a  one-layer  bounded  duct.  The  canonical  eigenvalue  equation  for 
a  multilayered  profile  is  a  function  of  two  dimensionless  mathematical  variables,  x  and  y,  and  two  sets 
of  dimensionless  parameters,  and  5/ .  The  parameter  p  j  is  a  function  of  the  ratio  of  sound-speed 
gradients  at  layer  interfaces,  while  the  parameter  5,  involves  ratios  of  three  successive  interface  sound 
speeds.  The  canonical  eigenvalues  are  curves  of  y  vs  x  for  a  fixed  set  of  p j  and  S,,  which  total  2(/-l) 
parameters  for  a  bounded  profile  of  /  layers. 

The  eigenvalues  are  referred  to  as  canonical  because  they  do  not  apply  to  just  one  profile  as  in 
the  customary  approach.  They  apply  for  some  given  boundaty  condition  to  an  entire  class  of  profiles, 
viz.,  the  subset  of  all  /-layered  profiles  as  constrained  by  the  fixed  set  of  p;-  and  5/. 

For  a  one-layered  bounded  profile,  there  are  no  p^  and  Si  parameters  to  constrain  the  profiles. 
The  set  of  curves  of  y  vs  x  for  a  given  boundary  condition  represents  the  canonical  eigenvalues  for  the 
entire  class  of  one-layer  bounded  profiles  with  the  given  boundary  condition.  This  report  deals  with 
three  pairs  of  boundary  conditions.  These  are  a  free  surface  and  rigid  bottom,  a  free  surface  and  free 
bottom,  and  a  rigid  surface  and  rigid  bottom. 

Section  1  presents  numerical  examples  for  canonical  eigenvalues,  phase  velocity,  group  velocity, 
and  eigenfunction  normalization  coefficients  as  well  as  for  certain  other  auxiliary  functions  referred  to 
as  canonical  frequency,  canonical  group  velocity  factor,  and  canonical  normalization  coefficient.  In 
the  computational  process,  these  auxiliary  functions  are  generated  along  with  the  canonical  eigen¬ 
values.  These  auxiliary  functions,  however,  will  not  be  presented  until  they  arise  in  the  orderly  pro¬ 
gression  of  results.  Numerical  examples  for  the  three  sets  of  boundary  conditions  are  presented  in  par¬ 
allel,  because  it  is  of  interest  to  compare  the  three  sets  as  the  progression  of  results  is  presented. 
Although  Section  1  is  primarily  concerned  with  numerical  results,  it  deals  with  theoretical  analysis  as 
the  need  arises. 

Section  2  presents  a  summary.  Appendixes  A,  B,  and  C  treat  the  more  involved  proofs  of  various 
characteristics  of  the  canonical  eigenvalues. 
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SECTION  1.  RESULTS  FOR  ONE-LAYER  BOUNDED  DUCTS 


The  one-layer  bounded  ducts  fall  into  three  profile  classes.  The  class  of  most  interest  in  underwa¬ 
ter  acoustics  is  illustrated  in  Fig.  1.  Here  the  surface  boundary  is  free,  while  the  bottom  boundary  is 
rigid.  As  with  all  canonical  profile  classes,  the  profiles  are  divided  into  two  cases:  Case  A  with  a  posi¬ 
tive  slope  and  Case  B  with  a  negative  slope. 

CASH  A  CASE  B 

FREE  SURFACE 


RIGID  BOTTOM 


Figuie  1.  Schematic  of  one-layer  duct  with  free  surface 
and  rigid  bottom. 

A  second  profile  class  is  that  of  the  free  surface,  free  bottom.  A  third  profile  class  is  that  of  the 
rigid  surface,  rigid  bottom.  These  two  profile  classes  are  primarily  of  interest  only  theoretically  for 
underwater  acoustics.  They  might,  however,  have  practical  application  in  other  disciplines  involving 
wave  equation  solutions.  Our  interest  here  is  how  these  different  boundary  conditions  affect  the  ca¬ 
nonical  eigenvalue  solutions  and  subsequent  results. 

/  fourth  profile  class— rigid  surface  and  free  bottom— presents  nothing  new,  since  the  results  are 
essentially  the  same  as  for  the  first  class  of  Fig.  1.  (If  we  turn  class  4  upside  down,  we  have  class  1.) 


We  next  define  four  functions  of  x  and  y: 

G,  =  Ai(-x)  Bi(-y)  -  Bi(-x)Ai(-y )  (1) 

G2  =  Ai{-x)  Bi'(-y)  -  Bi(-x)  Ai’(-y)  (2) 

G3  =  Ai'(-x)  Bi(-y )  -  Bi'(-x)  Ai(-y)  (3) 

G4  =  Ai’(-x)  Bi'(-y)  -  Bi'(-x)  Ai'(-y)  (4) 


Here  x  and  y  represent  the  Airy  function  argument  evaluated  at  the  surface  and  bottom,  respec¬ 
tively.  Most  of  the  functions  of  x  and  y  of  interest  in  the  canonical  approach  can  be  expressed  in 
terms  of  these  four  functions. 

CANONICAL  EIGENVALUES 

The  eigenvalue  equation  for  the  duct  of  Fig.  1  is 

Gi  =  0  (5) 

where  G2  is  given  by  Eq.  2.  Equation  5  guarantees  that  the  F i  of  Eq.  2  of  Ref.  1  vanishes  at  the  sur¬ 
face  boundary  and  that  dF\/dZ  vanishes  at  the  rigid  bottom.  Figure  2  presents  the  solutions  of  Eq.  5. 
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Observe  the  line  y  =  x.  As  discussed  in  Eq.  15  and  16  of  Ref.  1,  this  line  separates  the  modes  into 
two  families.  Those  modes  that  lie  below  the  line  are  associated  with  Case  A  of  Fig.  1;  the  modes  that 
lie  above  the  line  are  associated  with  Case  B.  As  we  shall  demonstrate  later,  the  numbers  associated 
with  Case  A  or  B  for  each  eigenvalue  curve  are  mode  numbers,  i.e.,  the  mode  count  increases  as 
each  family  member  moves  away  from  the  line  y  =  x. 


Figure  2  Canonical  eigenvalues  for  the 
free-rigid  boundaries. 


Consider  next  the  solid  vertical  and  horizontal  grid  lines.  The  vertical  lines  are  solutions  or  roots 
of 


Ai(-x)  =  0  (6) 

while  tl  *;zontal  lines  are  solutions  of 

Ai'(-y)  =  0  (7) 

The  numerical  value  of  these  roots  are  marked  on  the  grid  lines.  One  may  readily  verify  that  the  x,  y 
combinations  of  Eq.  6  and  7  satisfy  Eq.  5.  Thus,  the  intersections  of  the  grid  lines  represent  lattice 
points,  which  the  eigenvalue  curves  must  pass  through.  In  addition  to  the  finite  zeros,  Eq.  6  and  7  go 
to  zero  asymptotically  for  large  negative  values  of  x  and  y,  respectively.  The  lattice  points  at  x  =  -  °° 
are  manifest  in  Fig.  2  by  the  horizontal  asymptotes  at  the  zeros  of  Eq.  7  for  the  Case  B  curves.  The 
lattice  points  at  y  =  -  »  are  manifest  in  Fig.  2  by  the  vertical  asymptotes  at  the  zeros  of  Eq.  6  for  the 
Case  A  curves. 
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Now  consider  the  dashed  vertical  and  horizontal  grid  lines.  The  vertical  lines  are  solutions  of 


Bi(-x)  =  0  (8) 

while  the  horizontal  lines  are  solutions  of 

Bi'(-y)  =  0  W 

These  lines  are  also  labeled  with  the  numerical  value  of  the  roots.  These  x,  y  combinations  again  sat¬ 
isfy  Eq.  5  and  represent  lattice  points,  which  the  eigenvalue  curves  must  pass  through.  Equations  8 
and  9  have  only  finite  zeros.  Thus,  there  are  no  asymptotes  associated  with  the  dashed  grid  lines. 

One  can  readily  predict  the  curves  by  connecting  the  lattice  points  and  noting  that  each  lattice 
point  must  lie  on  one  curve.  The  only  areas  which  really  need  a  detailed  solution  of  Eq.  5  are  the 
case  A  curves  for  y  <  1.019  and  the  case  B  curves  for  x  <  1.174. 

The  curves  of  Fig.  2  were  generated  by  solving  Eq.  5  by  iteration  with  Newton’s  method.  Two 
approaches  were  used.  The  first  approach  was  to  select  a  y  and  iterate  on  x  to  find  the  x  solution  for 
the  fixed  y.  This  approach  makes  use  of 


6G2/6x  =  -G<  (10) 

The  second  approach  was  to  select  an  x  and  iterate  on  y  to  find  the  y  solution  for  the  fixed  x.  This 
apnroach  makes  use  of 


6G2/6y  =  yGi  (11) 

Our  approach  is  a  “mode  follower  routine”  which  traces  out  a  continuous  x,  y  curve  correspond¬ 
ing  to  some  mode.  We  start  the  iterative  process  near  a  lattice  point,  where  we  have  good  estimates  of 
the  solution,  and  trace  the  curve  by  moving  first  to  lower  values  of  x  and  y  and  then  to  higher  values 
of  x  and  y.  First  estimates  for  each  iteration  are  obtained  by  extrapolating  the  known  solution.  Once 
the  first  curve  has  been  traced,  we  can  trace  a  second  curve,  starting  at  any  lattice  point  that  does  not 
lie  on  the  first  curve.  Additional  curves  are  traced  by  starting  at  any  lattice  point  that  does  not  lie  on 
curves  already  traced. 

The  approach  of  Eq.  11  works  best  for  the  Case  B  curves;  that  of  Eq.  10  works  best  for  the  Case 
A  curves.  Equation  10  can  run  into  trouble  for  Case  B  curves,  because  the  tracing  process  may  select 
a  y  for  which  the  curve  does  not  exist.  Equation  1 1  can  run  into  trouble  for  Case  A  curves,  because 
the  tracing  process  may  select  an  x  for  which  the  curve  does  not  exist.  Moreover,  Eq.  11  vanishes  at 
y  -  0,  as  evidenced  by  the  vertical  slopes  for  the  Case  A  curves  at  this  point.  This  causes  the  routine 
to  fail  if  one  attempts  to  use  Eq.  1 1  in  the  region  near  y  =  0. 

We  have  already  associated  the  two  families  of  modes  with  Case  A  and  Case  B  of  Fig.  1  on  the 
basis  of  their  relation  to  the  line  y  =  x.  Ihere  is  a  second  mathematical  property  which  confirms  this 
association.  Observe  that  the  B  family  of  modes  all  cross  x  =  0,  but  not  y  =  0.  From  Eq.  17  of  Ref. 

1,  Cp  =  Ct  when  x  =  0.  This  corresponds  to  a  ray  which  grazes  the  ocean  surface  and  can  only 
occur  for  Case  B  in  Fig.  1.  Similarly,  the  A  family  of  modes  all  cross  y  =  0,  but  not  x  =  0.  From  Eq. 

8  of  Ref.  1,  we  can  establish  that  Cp  =  C2  when  y  =  0.  This  corresponds  ..o  a  ray,  which  grazes  the 
ocean  bottom,  and  can  only  occur  for  Case  A  in  Fig.  1. 
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We  next  address  the  free-surface,  free-bottom  profile.  Figure  3  is  the  counterpart  of  Fig.  2,  and 
the  counterparts  of  Eq.  5  to  1 1  are,  respectively 


C,  =  0 

(12) 

Ai{- x)  -  0 

(13) 

O 

II 

y— v 

1 

(14) 

Bi(-x)  =  0 

(15) 

Bi(-y)  =  0 

(16) 

6G\/dx  —  ~  G3 

07) 

6G\/dy  -  -G2 

(18) 

Figure  3.  Canonical  eigenvalues  for  the 
free-free  boundaries. 

Figure  2  has  many  features  similar  to  Fig.  1.  There  are,  however,  two  important  differences.  The 
first  difference  is  that  y  =  x  is  a  solution  to  Eq.  12.  This  is  readily  verified  from  Eq.  1.  Observe  also 
that  the  line  goes  through  a  set  of  lattice  points.  The  line  y  =  x  represents  a  degenerate  mathematical 
solution  which  is  discarded  physically  because  it  corresponds  to  zero  frequency  for  all  values  of  y  =  x. 
A  second  difference  is  that  the  two  families  of  modes  are  symmetric  about  the  line  y  =  x.  This  may 
be  demonstrated  by  noting  that  Eq.  12  is  unchanged  when  x  and  y  are  interchanged.  What  this  means 
physically  is  that  Case  B  may  be  considered  to  be  Case  A  turned  upside  down. 
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We  now  address  the  rigid-surface,  rigid-bottom  profile.  Figure  4  is  the  counterpart  of  Fig.  2,  and 
the  counterparts  of  Eq.  5  to  1 1  are,  respectively 


G4  =  0 

(19) 

Ai'(-x)  =  0 

(20) 

Ai'{-y )  =  0 

(21) 

Bi'{~  x)  =  0 

(22) 

Bi'(-y)  =  0 

(23) 

5G4/6x  -  xG2 

(24) 

6G4/dy  —  yGj 

(2V 

Figure  4.  Canonical  eigenvalues  for  the 
rigid-rigid  boundaries. 

Figure  4  is  similar  to  Fig.  3  in  that  the  line  y  =  x  is  a  solution  and  the  two  families  of  modes  are 
symmetric  about  the  line  y  =  x.  The  real  surprise  is  the  strange  behavior  of  mode  1.  Our  first  impres¬ 
sion  was  that  the  results  for  mode  1  were  in  error.  However,  as  discussed  later,  this  curve  is  compat¬ 
ible  with  the  result  that  mode  1  for  a  duct  with  rigid  boundaries  can  never  correspond  to  reflection 
from  both  boundaries.  In  Fig.  4,  the  mode  1  curves  meet  at  x  =  0.  y  =  0.  Mode  1A  never  crosses 
into  positive  y.  Thus,  the  mode  corresponds  to  a  ray  which  never  reflects  from  the  bottom  interface. 
Similarly,  Mode  IB  never  crosses  into  positive  x.  Thus,  the  mode  corresponds  to  a  ray  which  never 
reflects  from  the  surface  interface. 

The  nature  )f  the  mode  1  curve  can  be  inferred  from  »he  characteristics  of  the  lattice  points. 
There  must  be  a  IB  curve  which  is  asymptotic  to  y  =  1.019  and  a  1A  curve  which  is  asymptotic  to 
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x  =  1.019  because  of  the  zeros  of  Ai'  at  x  or  y  =  -  ».  However,  the  IB  and  1A  curves  cannot  cross 
into  y  >  1.019  or  t  ^  1.019  because  all  of  the  lattice  points  in  this  region  already  lie  either  on  one  of 
the  curves  for  higher-order  modes  or  on  the  line  y  =  x. 

W»  wul  need  to  know  the  nature  of  the  solution  for  1A  and  IB  as  the  origin  is  approached. 

Obsei  .e  in  Fig.  4  that  the  eigenvalue  curve  appears  to  cross  the  curve  y  =  x  at  right  angles.  This  sug¬ 
gests  that  the  line 


>’  =  -x  (26) 

is  an  approximate  solution  to  Eq.  19  for  mode  1.  If  we  assume  that  Eq.  26  holds,  C4  of  Eq.  4 
becomes 


G4( jc.-x)  =  Ai'(-x )  -  Bi'(-x)Ai'( x)  (27) 

We  next  evaluate  Eq.  27  by  means  of  the  series  expansions  given  in  Sections  10.4.2  to  10.4.5  of  Ref. 
2.  After  much  manipulation  and  simplification,  we  obtain 

Gt(x,-x )  =  (8/15)3  ,/2C,C2^5  (28) 

where  C\  and  C2  are  constants  given  in  Ref.  2.  The  numerical  value  of  the  coefficient  of  x5  in  Eq. 

28  is  0.0848826.  The  results  of  Eq.  27  were  compared  with  those  of  Eq.  28  for  x  =  0.5.  The  value  of 
Eq.  28  was  0.0026526  and  agreed  with  Eq.  27,  with  a  difference  of  1  in  the  fifth  significant  digit. 
Since  the  error  in  the  approximate  solution  of  Eq.  26  varies  as  x5 ,  it  is  accurate  near  the  origin. 
Indeed.  Fig.  4  shows  little  departure  from  Eq.  26  until  |  x  |  >  0.5.  Equation  26  will  be  useful  in  the 
analysis  of  phase  velocity  for  modes  1A  and  IB. 

There  is  another  feature  of  interest  in  Fig.  2  to  4.  Consider  the  spacing  between  modes  at  larger 
values  of  x  and  y,  away  from  the  vertical  and  horizontal  asymptotes.  In  Fig.  2,  the  spacing  is  about 
the  same  between  modes  1A  and  IB  as  it  is  between  adjacent  modes  of  family  A  and  family  B.  Simi¬ 
larly,  in  Fig.  3  and  4.  the  spacing  between  the  line  y  =  x  and  mode  1A  or  mode  IB  is  about  the 
same  as  it  is  between  adjacent  modes  of  family  A  and  family  B.  This  near-regular  spacing  between 
mathematical  solutions  is  a  property  of  solutions  to  eigenvalue  equations. 

Our  contention  is  that  the  eigenvalue  solutions  in  mathematical  space  are  more  predictable  and 
provide  more  information  than  in  customary  physical  space.  For  example,  the  mathematic  solutions 
for  both  mode  families  and  the  y  =  x  solutions  all  fit  into  a  consistent  pattern.  In  contrast,  the  cus¬ 
tomary  solutions  in  physical  space  for  family  A  have  no  apparent  connection  with  the  solutions  for 

family  B,  and  there  is  no  counterpart  at  all  for  the  solutions  given  by  y  =  x. 

Our  plans  call  for  a  systematic  investigation  of  canonical  eigenvalue  plots  for  two  and  even  more 

layers.  We  believe  that  these  plots  will  display  a  behavior  which  fits  into  a  regular  scheme  that  is  more 
predictable  than  the  customary  plots  in  physical  space. 
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CANONICAL  FREQUENCY  PLOTS 

Canonical  frequency  plots  apply  to  the  entire  class  of  profiles.  The  frequency  for  a  specific  profile 
of  the  class  is  given  by  Eq.  14  of  Ref.  1  and  may  be  obtained  by  multiplying  the  canonical  frequency 
by  a  scale  factor  which  depends  on  the  profile  parameters. 

Figures  5  to  7  are  the  respective  canonical  frequency  plot  counterparts  of  Fig.  2  to  4.  These  aie 
obtained  by  evaluating  Jjy  —  jc | 3/2  from  the  eigenvalue  curves  of  Fig.  2  to  4  and  plotting  the  result  vs  x. 
Mode  family  A  is  given  by  the  solid  curves  and  mode  family  B  by  the  dashed  curves.  These  curves 
apply  to  all  profiles  of  the  class. 

Before  discussing  the  propenies  o'  the  curves,  we  demonstrate  how  to  obtain  the  acoustic  fre¬ 
quency  for  some  particular  desired  profile  of  the  class.  We  shall  use  the  same  sound-speed  profiles 
throughout  this  report  when  we  need  to  evaluate  quantities  that  are  dependent  on  a  specific  profile. 
The  example  profile  chosen  for  Case  A  is  C\  =  1480  m/s,  Ci  -  1482  m/s,  and  y,  =  0.02s-1.  The 
example  profile  chosen  for  Case  B  is  C\  =  1482  m/s,  Ci  =  1480  m/s,  and  y,  =  -0.02s'1.  These  pro¬ 
files  are  selected  to  be  quite  similar  in  that  the  maximum  and  minimum  sound  speeds  are  the  same 
and  the  surface  gradient  has  the  same  magnitude.  The  layer  thickness  is  99.80  m  for  Case  A  and 
100.20  m  lor  Case  B.  The  frequency  scale  factor  of  Eq.  14  of  Ref.  1  is  45.44  for  Case  A  and  45.26 
for  Case  B.  Thus,  for  example,  the  frequency  corresponding  to  an  ordinate  of  4.0  in  Fig.  5  to  7  is 
about  182  Hz  for  the  given  profiles. 

Consider  now  the  properties  ot  the  canonical  frequency  curves.  As  x  approaches  °°,  the  frequency 
for  both  families  and  for  all  three  boundary  conditions  goes  to  zero  asymptotically.  A  proof  is  given  in 
Appendix  B.  where  it  is  demonstrated  tnat  y  -  x  goes  to  zero  for  all  modes.  What  happens  is  that  all 
of  the  eigenvalue  curves  of  Fig.  2  to  4  approach  the  line  y  =  x  asymptotically  for  large  values  of  x.  An 
exception  is  mode  1  cf  Fig.  4  and  7.  Here,  /  goes  to  zero  for  mode  1A  as  x  approaches  zero  from 
positive  values,  and  /  goes  to  zero  for  mode  IB  as  r  approaches  zero  from  negative  values.  This  zero 
limit  on  /  follows  immediately  from  Eq.  14  of  Ref.  1  and  Eq.  26. 

The  canonical  frequency  plots  always  have  vertical  asymptotes  for  the  A  family  of  modes  at  the 
same  values  of  *  as  the  vertical  asymptotes  of  the  corresponding  Fig.  2  to  4.  From  Eq.  14  of  Ref.  1, 
u.e  see  that  infinite  negative  y  and  fixed  x  leads  to  infinite  /.  For  the  Case  B  family,  there  are  no  ver¬ 
tical  asymptotes.  However,  from  Eq.  14  of  Ref.  1,  it  is  clear  that  /  goes  to  infinity  for  large  negative  x 
for  the  B  family  of  modes.  Here,  y  is  fixed  as  x  goes  to  minus  00 .  The  important  feature  here  is  that  / 
goes  to  infinity  for  either  x  -  »  or  y  .  Whether  the  curves  form  vertical  asymptotes  or  not  is 

not  important.  For  example  if  we  had  chosen  y  rather  than  x  as  the  independent  variable,  we  would 
find  that  the  B  family  formed  vertical  asymptotes  while  the  A  family  did  not. 
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Figure  5.  Canonical  frequency  plot  for 
the  free-rigid  boundaries. 
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Figure  6.  Canonical  frequency  plot  for 
the  free-free  boundaries. 


x 

Figure  7.  Canonical  frequency  plot  for 
the  rigid-rigid  boundaries. 
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PHASE  VELOCITY 


Figures  8  to  10  are  the  phase  velocity  counterparts  of  Fig.  2  to  4,  respectively.  These  curves  were 
generated  by  Eq.  17  of  Ref.  1  by  using  the  values  of  /  from  Eq.  14  of  Ref.  1,  the  canonical  eigen¬ 
value  x,  and  the  specified  values  of  C\  and  y\ .  Again,  family  A  values  are  given  as  solid  curves  and 
family  B  values  as  dashed  curves. 

We  note  that  the  A  and  B  curves  are  well  separated  in  Fig.  8.  However,  in  Fig.  9  and  10  the 
curves  appear  to  be  superimposed.  Appendix  A  demonstrates  that  if  the  phase  velocity  is  Cp  for  fre¬ 
quency  /  for  Case  A,  then  the  same  phase  velocity  occurs  for  frequency  /  for  Case  B,  where 

f  =  /(C,/C2)3  (29) 

Here,  Cx  and  C2  are  the  parameters  for  Case  A.  For  the  parameters  of  our  sample  profile, 

(C>/C2)3  =  0.9960.  Thus,  the  curves  for  Case  B  lie  very  slightly  to  the  left  of  those  for  Case  A. 

The  condition  of  Eq.  29  depends  on  the  symmetry  of  the  curves  of  Fig.  3  and  4  about  the  line 
y  =  x  and  also  in  our  choice  of  profile  parameters  for  Case  A  and  Case  B.  When  we  chose  the  pa¬ 
rameters  so  that  the  Case  A  and  Case  B  profiles  would  resemble  each  other,  we  had  no  idea  that  the 
phase  velocity  solutions  would  be  almost  identical  for  the  boundary  conditions  of  Fig.  3  and  4. 

Although  we  determined  phase  velocity  from  Eq.  17  of  Ref.  1,  the  following  alternative  expression 
is  useful  for  demonstrating  various  features  of  the  solution: 

Cp  =  C2[l  -  r1'3  *'2/3  |yio|2/3>]'1/2  (30) 

Equation  30  is  obtained  by  solving  Eq.  8  of  Ref.  1  for  Cp  with  /  =  1.  Both  Eq.  17  of  Ref.  1  and  Eq. 
30  hold.  This  follows  because  Eq.  14  of  Ref.  1  provides  the  condition  that  they  both  hold.  Indeed, 

Eq.  14  of  Ref.  1  is  most  easily  obtained  by  setting  Eq.  17  of  Ref.  1  and  Eq.  30  equal  and  solving  for 
/. 


Equation  30  makes  it  apparent  that  y  =  0  corresponds  to  Cp  =  C2,  as  was  previously  mentioned 
in  the  discussion  of  Fig.  2.  We  chose  to  consider  x  as  the  independent  and  y  as  the  dependent  vari¬ 
able  of  the  eigenvalue  equation.  However,  we  could  have  considered  y  as  the  independent  variable 
and  used  Eq.  30  to  determine  Cp  rather  than  Eq.  17  of  Ref.  1. 

Consider  now  the  characteristics  of  the  curves  of  Fig.  8  to  10  at  the  higher  frequencies.  From  Eq. 
17  of  Ref.  1,  it  is  clear  that  Cp  goes  to  Cx  for  the  A  family  of  modes,  since  x  is  a  finite  constant. 
Similarly,  Eq.  30  demonstrates  that  Cp  goes  to  C2  for  the  B  family  of  modes,  since  y  is  a  finite  con¬ 
stant.  In  Fig.  8  to  10,  at  high  frequencies,  Cp  goes  to  1480  m/s  for  both  families  of  modes,  because 
this  is  the  minimum  of  the  two  boundary  sound  speeds  chosen  for  our  sample  profiles. 

Observe  in  Fig.  8  to  10  that  the  phase  velocity  becomes  large  at  low  frequencies.  Indeed,  the 
phase  velocity  goes  to  infinity  for  specific  values  of  frequency.  These  frequencies  are  the  cutoff  fre¬ 
quencies  for  the  duct.  At  lower  frequencies,  the  phase  velocity  turns  pure  imaginary.  These  cutoff  fre¬ 
quencies  are  best  expressed  in  terms  of  wavelength  as  given  by 

/  =  No  Ci  Z2'  (31) 

where  N0  is  the  number  of  vertical  wavelengths  in  the  duct,  Ci  is  the  surface  sound  speed,  and  Z2  is 
the  duct  thickness. 


10 


Figure  8.  Example  of  phase  velocity  vs  Figure  9.  Example  of  phase  velocity  vs 

frequency  for  the  free-rigid  boundaries.  frequency  for  the  free-free  boundaries. 


Figure  10.  Example  of  phase  velocity  vs 
frequency  for  the  rigid-rigid  boundaries. 


The  number  of  wavelengths  at  cutoff  for  the  free-rigid  boundary  condition  is  given  by 


N0  =  (2n  -  l)/4  (32) 

where  n  is  the  mode  number.  The  corresponding  expression  for  the  free-free  and  rigid-rigid  boundary 
conditions  are,  respectively 


N0  =  n/2  (33) 

and 

N0  =  (n  -  l)/2  (34) 

A  rigorous  derivation  of  Eq.  31  to  34  is  given  in  Appendix  B.  The  derivation  is  laborious  and 
involves  asymptotic  representations  of  the  x,  y  lattice  points  and  several  applications  of  the  binomial 
expansion.  However,  the  results  can  be  explained  rather  simply  in  terms  of  the  trigonometric  func¬ 
tions,  which  represent  good  approximations  to  the  Airy  functions  for  these  arguments.  In  Eq.  32,  the 
standing  wave  is  like  a  sine  function,  with  a  node  at  the  surface  and  an  antinode  at  the  bottom.  Here, 
mode  n  has  n  antinodes  and  n  nodes.  The  antinodes  at  the  bottom  occur  at  odd  multiples  of  a  quar¬ 
ter  wavelength. 

In  Eq.  33,  the  standing  wave  is  like  a  sine  function  with  nodes  at  the  surface  and  bottom.  Here, 
mode  n  has  n  antinodes  and  n  +  1  nodes.  The  nodes  at  the  bottom  occur  at  multiples  of  one-half 
wavelength.  In  Eq.  34,  the  standing  wave  is  like  a  cosine  function,  with  antinodes  at  the  surface  and 
bottom.  Here,  mode  n  has  n  antinodes  and  n  -  1  nodes.  The  antinodes  at  the  bottom  occur  at  multi¬ 
ples  of  one-half  wavelength. 

We  note  that  for  the  sample  profile  parameters 

CJZi  =  14.83  or  14.79.T1  (35) 

or  Case  A  or  Case  B,  respectively.  Thus,  the  cutoff  frequencies  for  Case  A  or  B  are  at  odd  multiples 
of  3.71  or  3.70  Hz,  respectively,  in  Fig.  8  and  are  at  all  multiples  of  7.42  or  7.40  Hz  in  Fig.  9  and 
10. 


On  our  first  exposure  to  the  cutoff  frequencies,  the  result  appeared  at  odds  with  a  fundamental 
property  for  Sturm-Liouville  equations,  as  discussed  in  Ref.  3.  For  our  boundary  conditions,  the 
Sturm-Liouville  eigenvalues,  k,  are  real;  whereas  our  eigenvalues,  Cp,  turn  imaginary  below  cutoff. 
However,  A  =  a)2/Cp .  Thus,  since  our  values  of  Cp  are  not  complex  but  pure  imaginary  below  cutoff, 
there  is  no  contradiction.  The  eigenvalue  X.  is  positive  for  frequencies  above  cutoff  and  negative  for 
frequencies  below  cutoff. 

We  turn  now  to  a  discussion  of  mode  1  for  the  rigid-rigid  boundary.  There  is  no  mode  1  for  sur¬ 
face-reflected  rays.  This  was  pointed  out  in  the  discussion  of  Fig.  4  and  is  in  agreement  with  the 
results  of  the  phase-integral  method  as  discussed  in  Ref.  4.  For  mode  1,  Eq.  40  of  Ref.  4  is  zero  and 
the  phase-integral  method  has  no  solution.  Also  in  Ref.  4,  the  phase-integral  solution  for  mode  n  of 
the  free-free  boundary  is  the  same  as  for  mode  n  -  1  of  the  rigid-rigid  boundary.  We  may  determine 
from  Fig.  3  and  4  that  this  is  a  very  good  approximation.  The  lattice  points  for  n  in  Fig.  3  are  similar 
to  those  for  mode  n  +  1  in  Fig.  4.  This  is  true  because  the  zeros  of  Ai  lie  close  to  those  of  Bi'  and 
the  zeros  of  Bi  lie  close  to  those  of  Ai' . 

For  an  isovelocity  duct,  there  is  no  mode  1  for  the  rigid-rigid  boundary  condition.  However,  for 
the  profile  of  Eq.  1  of  Ref.  1  with  a  sound-speed  gradient,  there  is  a  mode  1  solution  for  refracted 


12 


rays,  which  reflect  from  only  one  boundary.  The  phase-integral  method  gives  an  approximate  solution 
for  n  =  1  in  Eq.  39  of  Ref.  4.  We  next  examine  the  behavior  of  phase  velocity  for  mode  1  as  i  goes 
to  zero.  Equation  26  holds  here,  and  the  /  of  Eq.  14  of  Ref.  1  goes  to  zero.  The  substitution  of  Eq. 
26  into  Eq.  14  of  Ref.  1  and  hence  into  Eq.  17  of  Ref.  1  leads  to 

Cp  =  [(Cl2  +  c;2)/  2]->/2  (36) 

In  other  words,  the  square  of  the  phase  velocity  corresponding  to  /  =  0  for  mode  1  of  the  rigid-rigid 
boundary  condition  is  the  harmonic  mean  of  the  square  of  the  two  boundary  sound  speeds. 

In  Fig.  10,  the  phase  velocity  for  /  =  0  appears  to  be  1481— midway  between  the  boundary  sound 
speeds.  If  we  let 


C2  =  C,  +A 


(37) 


we  may  expand  Eq.  36  as  a  series  in  A  to  yield 

Cp  s  C,  +  A/2  -  3A2/8C,  (38) 

The  value  of  the  first  two  terms  of  Eq.  38  is  midway  between  Cx  and  C2.  For  the  sample  profiles,  the 
third  term  of  Eq.  38  reduces  Cp  from  the  midpoint  value  to  1480.9990  m/s. 

We  note  that  Eq.  34  gives  the  correct  result  for  n  =  1  in  that  the  cutoff  frequency  for  mode  1  is 
zero.  However,  this  result  appears  fortuitous  because  Eq.  34  represents  the  frequencies  where  the 
phase  velocity  is  infinite,  while  /  equals  zero  results  in  the  finite  value  of  Eq.  36.  Not  only  is  there  no 
mode  1  for  surface-reflected  rays;  there  is  no  mode  1  for  refracted  rays  with  phase  velocities  between 
Eq.  36  and  the  larger  of  the  sound  speeds  at  the  two  boundaries. 


GROUP  VELOCITY 

Figures  11  to  13  present  the  canonical  group  velocity  factor  L  of  Eq.  67  of  Ref.  1,  corresponding, 
respectively,  to  the  eigenvalue  curves  of  Fig.  2  to  4.  The  discussion  of  Fig.  11  to  13  is  best  included 
as  part  of  the  discussion  of  the  group-velocity-vs-frequency  curves,  which  are  shown  in  Fig.  14  to  16, 
respectively. 

These  curves  are  generated  from  Eq.  69  of  Ref.  1,  which  is  interesting  in  that  the  frequency  is 
not  contained  explicitly,  but  is  implicitly  in  Cp  and  L.  The  procedure  is  to  evaluate  L  and  /  from  the 
eigenvalue  curves,  then  evaluate  Cp  for  that  frequency,  and  then  evaluate  Cg  for  the  Cp  and  L  that 
corresponds  to  the  frequency. 

As  was  the  case  with  the  phase  velocity,  the  group-velocity  curves  of  Case  A  and  Case  B  are  al¬ 
most  superimposed  in  Fig.  15  and  16.  With  the  use  of  Eq.  63  of  Ref.  1  and  Eq.  29,  we  can  readily 
show  that  the  group  velocity  is  the  same  function  of  phase  velocity  for  Case  A  as  for  Case  B.  This 
property  can  also  be  demonstrated  from  Eq.  69  of  Ref.  1,  but  involves  a  tangle  of  involved  algebraic 
manipulation  which  incredibly  simplifies  to  the  result  obtained  from  Eq.  63  of  Ref.  1  and  Eq.  29.  The 
net  result  is  similar  to  that  for  the  phase  velocity.  If  the  group  velocity  is  Cg  for  frequency  /  for  Case 
A,  the  same  group  velocity  occurs  for  frequency  /  for  Case  B,  where  /  is  given  by  Eq.  29.  Again, 
the  group-velocity  curves  for  Case  B  lie  very  slightly  to  the  left  of  those  for  Case  A. 


13 


Figure  11.  Canonical  group  velocity 
factor  for  the  free-rigid  boundaries. 
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Figure  12.  Canonical  group  velocity 
factor  for  the  free-free  boundaries. 
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Figure  15.  Example  of  group  velocity  vs 
frequency  for  the  free-free  boundaries. 
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Figure  16.  Example  of  group  velocity  vs 
frequency  for  the  rigid-rigid  boundaries. 


It  seems  somewhat  incredible  that  Eq.  69  of  Ref.  1  can  transform  the  obscure  mathematical 
curves  of  Fig.  11  to  13  into  the  familiar  physical  curves  of  Fig.  14  to  16.  However,  we  next  provide  a 
physical  interpretation  of  Eq.  69  of  Ref.  1.  We  consider  first  the  Case  A  configuration.  As  indicated 
in  the  previous  discussion  of  Eq.  69  and  70  of  Ref.  1,  the  solution  for  the  unbounded  ducts  of  Ref.  4 
is  Eq.  69,  with  the  third  term  in  the  brackets  suppressed.  This  third  term  may  then  be  considered  as 
a  correction  term  to  the  group  velocity  as  given  by  the  familiar  Eq.  70  of  Ref.  1. 

Consider  now  the  Case  A  ducts  at  the  vertical  asymptotes  of  Fig.  2  to  4.  The  slope  dy/dx  is  +>» 
for  Fig.  3  and  -®  for  Fig.  2  and  4.  Thus,  L  of  Fig.  12  starts  at  -0  at  the  x  of  the  asymptote  and  at 
+0  in  Fig.  11  and  13.  For  Case  A,  the  sign  of  the  correction  term  in  Eq.  69  of  Ref.  1  is  the  negative 
of  that  ot  L.  Consider  now  the  free  bottom  of  Fig.  12.  Here,  L  starts  out  at  -0  and  continues  to  grow 
large  negatively.  Thus,  the  third  (correction)  term  starts  out  at  zero  and  increases  monotonically.  This 
in  turn  causes  the  group  velocity  to  be  smaller  than  the  result  of  Eq.  70  of  Ref.  1. 

The  behavior  is  a  bit  more  complicated  for  the  rigid  bottom  of  Fig.  11  and  13.  Here,  L  starts  at 
+0  and  forms  a  positive  relative  maximum  at  the  points  of  inflection  of  the  curves  of  Fig.  2  and  4.  It 
then  crosses  zero  at  y  =  0.  Thus,  for  the  rigid  bottom,  the  correction  term  is  slightly  negative  for 
refracted  rays,  i.e.,  y  <  0,  and  then  turns  increasingly  positive  for  surface-reflected  rays.  The  result  is 
that  the  group  velocity  for  refracted  rays  is  generally  somewhat  larger  than  Eq.  70  of  Ref.  1  and  is 
less  than  Eq.  70  for  surface-reflected  rays.  Here,  there  are  two  cases  where  Eq.  70  of  Ref.  1  gives  the 
group  verity,  at  infinite  irequency  corresponding  to  the  vertical  asymptote  and  for  the  ray  which 
grazes  the  rigid  bottom. 

The  foregoing  analysis  does  not  apply  to  mode  1A  for  the  rigid-rigid  condition  of  Fig.  13.  Here, 
dy/dx  goes  to  -1  at  x  =  0,  y  =  0  and  L  goes  to  0.5.  In  Fig.  10,  we  see  that  the  phase-velocity  curve 
has  zero  slope  at  zero  frequency.  Thus,  from  Eq.  63  of  Ref.  1,  we  determine  that  in  this  case 

Cg  =  Cp  (39) 

or  the  group  velocity  in  Fig.  16  is  given  by  Eq.  36.  One  may  also  obtain  Eq.  39  from  Eq.  69  of  Ref. 

1  by  noting  that  L  =  1/2  and  that  Eq.  36  holds. 

Consider  now  the  Case  B  modes.  Here,  we  can  express  Eq.  69  of  Ref.  1  as 

Cg  =  Cp[  \  +  (Cp  -  C|)/3Cl  +  C\(C ?  -  C\Z)(L  -  l)/3]-‘  (40) 

For  the  Case  B  profile  configuration,  one  must  consider  a  negative  rather  than  positive  half  space, 
and  the  roles  of  Cx  and  C2  are  reversed.  For  the  Case  B  modes,  Eq.  70  of  Ref.  1  corresponds  to  Eq. 
40,  with  the  third  term  in  the  brackets  suppressed.  The  third  term  of  Eq.  40  may  now  be  regarded  as 
a  corrective  term  to  the  result  of  Eq.  70  of  Ref.  1. 

At  the  horizontal  asymptotes  of  Fig.  2  to  4 ,  L  approaches  1  from  above  in  Fig.  2  and  3  for  the 
free  surface  and  from  below  in  Fig.  4  for  the  rigid  surface.  For  the  free  surface,  L  -  1  monotonically 
increases  and  the  group  velocity  is  smaller  than  Eq.  70  of  Ref  1.  This  is  the  same  behavior  as  for 
Case  A  for  the  free  bottom.  The  group  velocity  is  given  by  Eq.  70  of  Ref.  1  for  infinite  frequency, 
but  falls  progressively  below  the  result  of  Eq.  70  as  the  frequency  is  decreased. 

For  the  rigid  surface  of  Fig.  4  and  13,  L  -  1  is  negative  for  x  <  0.  There  is  a  relative  minimum  in 
L  -  1  at  the  point  of  inflection  in  Fig.  4,  and  it  again  crosses  zero  at  x  =  0  and  continues  to  grow 
positive  without  bound  as  x  is  increased.  Here  again,  the  behavior  is  the  same  as  for  Case  A  for  the 
rigid  bottom.  The  group  velocity  equals  that  ot  Eq.  70  of  Ref.  1  for  infinite  frequency  and  for  the 
frequency  corresponding  to  the  ray  which  grazes  the  rigid  bottom.  Between  these  frequencies,  the 
group  velocity  is  larger  than  that  of  Eq.  70  of  Ref.  1.  For  frequencies  corresponding  to  surface- 
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reflected  rays,  the  group  velocity  falls  progressively  below  Eq.  (70)  of  Ref.  1  as  the  frequency  is 
decreased. 

Again,  the  foregoing  analysis  does  not  apply  to  mode  IB  for  the  rigid-rigid  condition  of  Fig.  16. 
Here,  L  does  not  form  a  relative  minimum,  but  decreases  monotonically  with  increasing  x  to  join  with 
the  mode  1A  result  at  x  =  0  and  y  =  0.  Here,  Eq.  39  applies  for  zero  frequency. 

Consider  now  some  of  the  characteristics  of  Fig.  14  to  16.  At  high  frequencies,  Cg  asymptotically 
approaches  Cp.  This  follows  from  Eq.  69  of  Ref.  1  for  Case  A  modes  as  Cp  -*  C\  and  L  -*  0-  It  fol¬ 
lows  from  Eq.  40  for  Case  B  modes  as  Cp  -*  C\  and  L  -*•  1 .  The  limit  of  Cg  for  both  Case  A  and  B 
modes  is  1480  m/s  and  follows  from  our  choice  of  sample  profiles. 

As  the  frequency  is  decreased  from  large  values,  the  group  velocity  increases  and  forms  a  relative 
maximum,  the  so-called  Airy  phase  phenomenon,  and  then  decreases  rapidly  as  the  frequency  is  low¬ 
ered  further.  At  the  duct-cutoff  frequencies,  the  group  velocity  goes  to  zero.  This  can  be  verified  from 
Eq.  69  of  Ref.  1  as  Cp  -*  «  at  the  cutoff  frequencies,  and  the  third  term  in  the  brackets  is  positive. 

Again,  mode  1A  in  Fig.  16  has  a  different  behavior.  Here,  there  is  a  monotonic  increase  with 
decreasing  frequency,  with  an  end-point  maximum  at  zero  frequency  given  by  Eq.  39  and  36. 

The  effect  of  the  boundary  condition  on  the  group  velocity  is  best  illustrated  by  Fig.  14,  where  we 
see  that  the  group  velocity  for  Case  B  modes  with  the  free  surface  lies  well  below  the  group  velocity 
for  Case  A  modes  with  the  rigid  surface.  We  can  attribute  this  difference  to  two  effects.  We  note  that 
for  values  of  Cp  that  are  comparable  to  C\  and  C2>  Cp  appears  as  a  multiplicative  factor  in  Eq.  69  of 
Ref.  1.  Thus,  the  first  effect  is  the  lower  phase  velocities  in  Fig.  8  for  the  Case  B  modes.  We  have 
already  pointed  out  that  the  group  velocity  is  greater  than  Eq.  70  of  Ref.  1  for  refracted  rays  for  the 
rigid  boundary.  However,  it  is  smaller  than  Eq.  70  of  Ref.  1  for  refracted  rays  for  the  free  bounda¬ 
ries.  This,  then,  is  the  second  effect  producing  lower  group  velocities  for  the  free  boundaries. 

The  tendency  for  the  free  boundary  to  lower  group  velocities  is  also  evident  in  Fig.  15,  where  the 
group  velocity  for  both  sets  of  modes  is  a  relatively  flat  function  of  frequency,  with  the  maximum 
group  velocity  rising  slightly  above  1480  m/s.  In  contrast,  the  group  velocity  for  both  sets  of  modes  in 
Fig.  16  with  rigid-rigid  boundaries  peaks  significantly  higher  above  1480  m/s. 

EIGENFUNCTION  NORMALIZATION  CONSTANTS 

Before  presenting  numerical  results,  we  need  to  simplify  the  general  formulation  of  Ref.  1  so  it 
will  apply  to  a  single  layer.  We  will  now  present  the  theoretical  expressions  for  all  boundary  condi¬ 
tions.  For  a  single  layer  profile,  Eq.  90  of  Ref.  1  reduces  to 

Dc  =  Ds  +  Db  (41) 

Consider  the  free-rigid  duct  of  Fig.  1.  With  the  use  of  Eq.  81  of  Ref.  1,  the  F,(Z)  of  Eq.  2  of 
Ref.  1  may  be  written  as 

F,(Z)  =  Bi(-x)  Ai(-x  -  a,Z)  -  Ai(-x)  Bi(-x  -  a,Z)  (42) 

From  Eq.  73  and  81  of  Ref.  1  and  Eq.  1,  we  obtain 

F«(-f  io)  =  F,(-y)  =  -G,  (43) 

Equations  94  and  99  of  Ref.  1  and  Eq.  41  and  43  lead  to 

Dc  =  -Jr'2  t  y(-G|)2  (44) 
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Consider  next  the  free-free  duct.  Equation  42  applies  as  is.  From  Eq.  73  and  76  of  Ref.  1  and 
Eq.  2,  we  obtain 

to)  =  Ft'(-y)  =  G2  (45) 

Equations  94  and  98  of  Ref.  1  and  Eq.  41  and  45  lead  to 

Dc  =  -Jr'2  +  G\  (46) 

Consider  now  the  rigid-rigid  duct.  With  the  use  of  Eq.  82  of  Ref.  1,  the  F,(Z)  may  be  written  as 
Fi(Z)  =  Bi'(-x)  Ai{-x  -  a,Z)  -  Ai'(-x)  Bi(-x  -  a{Z)  (47) 

From  Eq.  73  and  82  of  Ref.  1  and  Eq.  3,  we  obtain 

F,(~U)  =  Fi(-y)  =  ~G3  (48) 

Equations  95  and  99  of  Ref.  1  and  Eq.  41  and  48  lead  to 

Dc  =  -xn'2  +  y  (~G3)2  (49) 

Figures  17  to  19  present  the  respective  results  of  Eq.  44,  46,  and  49  for  the  canonical  normaliza¬ 
tion  coefficient.  We  will  be  pointing  out  various  limits  as  x  or  y  becomes  large,  whether  negative  or 
positive.  The  mathematical  demonstrations  of  these  limits  are  quite  involved  and  are  presented  in 
Appendix  C. 

Note  first  that  with  one  exception  the  curves  never  cross  zero.  This  result  was  anticipated,  because 
Dc  equal  to  zero  would  correspond  to  the  case  where  a  mode  is  orthogonal  to  itself.  This  could 
happen  in  the  case  of  degenerate  eigenvalues,  which  cannot  occur  for  a  single  layer  duct.  In  point  of 
fact,  we  were  quite  relieved  when  Fig.  17  to  19  were  first  generated,  because  had  Dc  crossed  zero, 
we  would  have  known  that  we  had  made  an  error  in  developing  the  theoretical  expressions.  Note,  too, 
that  the  values  for  the  A  and  B  curves  are  always  negative  and  positive,  respectively. 

The  exceptional  curves  are  1A  and  IB  in  Fig.  19,  which  hook  together  at  Dc  =0  and  x  =  0. 

This  corresponds  to  zero  frequency.  However,  as  will  be  discussed  presently,  this  result  leads  to  a 
nonzero  value  for  Dn  and  apparently  poses  no  problem. 

All  of  the  A  and  B  curves  approach  zero  asymptotically  for  large  positive  x.  For  large  negative  x, 
the  B  curves  appear  to  increase  without  bound.  Note  that  the  A  curves  end  abruptly  at  finite  values  of 
Dc  and  x.  The  horizontal  line  in  Fig.  17  to  19  is  -  n~l-  The  x  corresponds  to  the  vertical  asymptotes 
in  Fig.  2  to  4.  The  G\,  G2,  and  G3  of  Eq.  44,  45,  and  49,  respectively,  are  messy  indeterminant 
forms  which  go  to  zero  in  the  limit.  As  a  result,  the  limit  of  Dc  goes  to  -n'1  in  Fig.  17  and  18,  and 
to  - n ~2  times  the  zeros  of  Ai'  in  Fig.  19.  Mode  1A  of  Fig.  19  does  not  end  on  the  horizontal  line  as 
do  all  the  A  modes  of  Fig.  17  and  18,  but  it  comes  close,  because  the  first  zero  of  Ai'  is  equal  to 
1.018.  In  the  initial  computer  runs,  these  finite  limits  on  the  A  curves  were  somewhat  ragged.  Double¬ 
precision  calculations  were  then  implemented  to  produce  the  clear-cut  results  of  Fig.  17  to  19. 

Figures  20  to  22  present  the  eigenvalue  normalization  coefficient  for  the  respective  curves  of  Fig. 
17  to  19.  These  curves  were  generated  by  Eq.  89  of  Ref.  1  by  using  the  specific  profile  parameters 
already  discussed.  Consider  first  the  A  curves,  which  appear  well-behaved.  Appendix  C  demonstrates 
that  they  all  go  to  zero  for  infinite  frequency.  Note  that  the  curves  of  Fig.  20  and  21  appear  to  merge 
to  a  common  curve  for  high  frequency.  This  is  an  intriguing  result,  but  we  are  not  inclined  to  tackle 
the  theoretical  analysis  in  this  article. 
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Figure  18.  Canonical  normalization 
coefficient  for  the  free-free  boundaries. 
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Figure  22.  Eigenvalue  normalization 
coefficient  for  the  rigid-rigid  boundaries. 


Consider  now  the  limits  of  D„  as  the  frequency  is  decreased  to  small  values.  Appendix  C  demon¬ 
strates  that  for  the  free-rigid  and  free-free  boundaries  of  Fig.  20  and  21,  respectively,  Dn  goes  to 
zero  at  the  duct  cutoff  frequencies.  However,  for  the  rigid-rigid  boundary  of  Fig.  22,  Dn  for  modes 
1A  and  IB  at  zero  frequency  goes  to  a  finite,  nonzero,  value  of  Z2;r"2.  where  Z2  is  the  duct  thick¬ 
ness.  This  result  is  also  demonstrated  in  Appendix  C.  The  limit  of  Dn  at  the  cutoff  frequencies  for 
both  A  and  B  curves  for  mode  2  and  higher-order  modes  is  apparently  Z2n'2/1-  We  say  apparently 
because  a  rigorous  proof  of  this  result  is  extremely  difficult,  and  the  discussion  of  this  result  in  Appen 
dix  C  is  theoretically  on  somewhat  shaky  ground. 

The  numerical  values  of  the  examples  of  Fig.  22  are  5.06  and  5.08  m  for  modes  2  and  above  for 
Cases  A  and  B,  respectively,  with  double  these  values  for  modes  1A  and  IB,  respectively.  Difficulties 
with  the  analysis  of  Appendix  C  led  us  to  question  at  one  time  whether  the  numerical  results  for 
modes  2A,  3A,  and  4A  in  Fig.  22  were  correct.  However,  the  result  for  1A  and  IB  is  on  firm  ground 
theoretically,  and  since  the  relative  maximum  of  2 A,  3 A,  and  4 A  line  up  well  with  the  1A  result,  we 
cannot  fault  the  numerical  result.  The  results  point  up  the  fact  that  Dn  in  Fig.  20  to  22  is  a  quantity 
with  dimensions  of  meters.  If  one  carries  out  the  process  forward  to  the  point  of  evaluating  propaga¬ 
tion  losses,  these  losses  are  then  referenced  relative  to  1  m. 

Whereas  the  normalization  process  appears  to  work  well  for  the  A  family  of  modes,  the  rapid 
increase  with  frequency  for  the  B  family  of  modes  makes  this  approach  unattractive.  We  were 
prompted  to  examine  an  alternative  method  of  normalization  based  on  the  boundary  condition  at  the 
bottom  rather  than  at  the  surface.  An  alternative  expression  for  Eq.  42  is 

FdZ)  =  Bi'(-y)  Ai(-x  +  a,Z)  -  Ai'(-y)  Bi(-x  +  a,Z)  (50) 

Here,  the  coefficients  Dx  and  £,  have  been  chosen  to  satisfy  the  bottom  boundary  condition  for  a 
rigid  bottom.  Equation  50  also  applies  to  the  rigid-rigid  duct.  The  counterpart  of  Eq.  50  for  the  free- 
free  duct  is 

Fi(Z)  =  Bi(-y)  Ai (- x  +  axZ)  -  Ai(-y)  Bi(-x  +  a,Z)  (51) 

With  this  alternative  method,  the  counterparts  of  Eq.  44,  46,  and  49  become,  respectively 

Dc  =  yn1  -  (-G4)2  (52) 

Dc  =  * -2  -  (-G3)2  (53) 

and 

Dc  =  yn1  -  x  G l  (54) 

Figures  23  to  25  present  the  respective  results  of  Eq.  52  to  54.  In  common  with  Fig.  17  to  19, 
the  values  of  Dc  are  positive  for  the  B  family  and  negative  for  the  A  family.  The  results  also  go  to 
zero  for  large  values  of  x.  In  contrast  to  Fig.  17  to  19,  the  behavior  of  the  A  and  B  families  is  re¬ 
versed.  Here,  the  B  family  goes  to  fixed  limits  for  large  negative  x,  while  the  A  family  diverges  by 
forming  vertical  asymptotes  as  y  goes  to  minus  infinity  in  Fig.  2  to  4.  Here,  the  G4 ,  G3 ,  and  G2  of 
Eq.  52  to  54,  respectively,  go  to  zero  in  the  limit.  As  a  result,  the  limit  of  Dc  goes  to  nl  times  the 
zeros  of  Ai'  in  Fig.  23  and  25  and  to  n~2  in  Fig.  24. 


23 


-4  -2  0  2  4  6  8 

X 

Figure  23.  Counterpart  of  Fig.  17 
for  alternative  normalization. 


Figure  24.  Counterpart  of  Fig.  18 
for  alternative  normalization. 
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Figure  25.  Counterpart  of  Fig.  19 
for  alternative  normalization. 


Figures  26  to  28  present  the  eigenvalue  normalization  coefficient  for  the  respective  curves  of  Fig. 
23  to  25.  We  see  that  the  alternative  approach  to  normalization  has  achieved  the  desired  goal.  Here, 
the  Case  B  modes  are  all  well  behaved,  although  the  Case  A  modes  now  increase  rapidly  with  increas¬ 
ing  frequency.  The  opposite  behavior  of  Case  A  and  B  modes  is  discussed  further  in  Appendix  C. 

In  closing  this  section,  we  note  that  the  dashed  and  solid  curves  of  Fig.  21  are  almost  identical  to 
the  solid  and  dashed  curves,  respectively,  of  Fig.  27.  The  same  property  holds  for  Fig.  22  and  28. 

The  choice  of  profile  parameters  is  such  that  if  we  turn  the  free-free  or  rigid-rigid  profiles  upside 
down,  we  have  almost  the  same  physical  problem,  with  the  roles  of  the  A  and  B  families  reversed. 

The  same  physical  problem  should  produce  the  same  results,  despite  the  fact  that  the  two  approaches 
are  quite  different.  We  will  not  attempt  to  demonstrate  theoretically  that  the  two  formulations  lead  to 
the  same  result,  because  such  a  demonstration  appears  very  involved.  For  example,  there  is  small 
resemblance  between  the  intermediate  results  of  Fig.  19  and  25.  Yet  the  final  results  of  Fig.  22  and 
28  agree.  We  were  pleased  at  the  result,  because  it  tends  to  increase  our  confidence  that  the  normali¬ 
zation  process  has  been  properly  implemented.  For  the  phase-  and  group-velocity  plots  of  Fig.  8  to  10 
and  Fig.  14  to  16,  we  had  enough  experience  with  the  customary  eigenvalue  approach  to  recognize 
that  the  canonical  approach  was  producing  familiar  patterns.  In  contrast,  we  had  no  such  experience 
with  plots  of  Dn  as  a  function  of  frequency,  and  hence  had  little  feel  for  what  the  characteristics  of 
Fig.  20  to  22  or  26  to  28  should  be. 
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Figure  28.  Counterpart  of  Fig.  22  for 
alternative  normalization. 


As  this  report  was  about  to  be  submitted  for  publication,  some  numerical  results  on  two-layer 
bounded  profiles  became  available.  These  results  suggested  that  the  normalization  curves  for  the  B 
modes  in  Fig.  17  to  19  may  not  diverge  for  negative  values  of  x.  The  case  of  Fig.  17  was  investigated 
in  more  detail,  and  it  was  found  that  for  mode  1  a  relative  maximum  in  Dc  occurred  at  x  =  -8.9, 
with  a  value  of  Dc  =  3  8  x  1014.  Furthermore,  this  result  implies  that  Dn  in  Fig.  20  to  22  does  not 
diverge  at  large  frequencies  but  forms  a  relative  maximum.  The  result  indicates  that  the  limit  analysis 
attempted  in  Appendix  C  is  even  more  difficult  than  believed  when  Appendix  C  was  prepared.  The 
question  of  whether  Dc  and  D„  diverge  for  the  B  modes  is  somewhat  moot,  since  the  alternative  nor¬ 
malization  of  Eq.  50  to  54  eliminates  the  problem. 
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SUMMARY 


The  application  of  the  theory  of  canonical  eigenvalues  to  underwater  sound  propagation  in  a  one- 
layer  bounded  duct  has  more  than  lived  up  to  the  expectations  we  had  for  it.  This  application  demon¬ 
strates  not  only  the  general  advantages  cited  in  Ref.  1,  but  also  some  new  advantages.  It  has  also  pro¬ 
duced  a  number  of  specific  results  of  importance  or  interest.  This  section  is  a  summary  of  these 
advantages  and  specific  results,  presented  for  the  most  part  in  the  order  in  which  they  arose  in  the 
approach. 

The  method  reduces  the  number  of  parameters  in  the  customary  eigenvalue  equation  by  four.  The 
four  parameters  eliminated  are  Ct,  C2<  Y\  anc*  /■  For  a  one-layer  bounded  duct,  this  leaves  no 
parameters  in  the  canonical  eigenvalue  equation.  Thus,  the  eigenvalue  solutions  apply  to  all  possible 
sound-speed  profiles  with  the  given  boundary  condition. 

The  canonical  eigenvalue  curves  pass  through  two  sets  of  lattice  points,  which  are  known  solutions 
to  the  eigenvalue  equation  and  are  obtained  from  zeros  to  various  Airy  functions.  For  reflection  from 
both  layer  boundaries,  these  lattice  points  almost  delineate  the  eigenvalue  curves.  Of  particular  impor¬ 
tance  are  the  zeros  at  minus  inanity,  which  define  the  vertical  asymptotes  for  the  Case  A  modes  and 
the  horizontal  asymptotes  for  the  Case  B  modes.  These  asymptotes  provide  important  characteristics  of 
the  propagation  at  high  frequencies  for  refracted  rays  which  reflect  from  only  one  boundary. 

The  near-uniform  spacing  between  the  set  of  A  and  B  modes  (and  of  the  line  y  -  x  when  it  is  a 
solution)  is  an  example  of  the  mathematical  essence  of  the  solution.  This  spacing  provides  a  clue  to 
the  location  of  mode  1,  which  is  not  available  in  the  customary  approach. 

The  plots  of  phase  velocity  vs  frequency  exhibited  the  familiar  characteristic  of  diverging  at  low 
frequency  and  asymptotically  approaching  the  lower  of  the  boundary  sound  speeds  at  high  frequency. 
This  divergence  is  related  to  a  low-frequency  cutoff,  which  detailed  analysis  demonstrated  could  be 
given  by  duct  thickness,  expressed  in  terms  of  acoustic  wavelength. 

The  expression  for  group  velocity  contains  the  sum  of  three  terms.  The  first  two  terms  can  be 
associated  with  the  classic  result  of  ray  theory,  while  the  third  term  is  a  “tweaking"  term  which  is  a 
function  of  the  slope  of  the  eigenvalue  curves.  This  expression  is  not  only  appealing  from  a  theoretical 
standpoint,  it  is  also  simpler  to  apply  than  the  formulation  used  in  the  customary  normal-mode 
approach.  The  group  velocities  for  refracted  rays,  forming  a  vertex  below  a  free  surface  or  a  nadir 
above  a  free  bottom,  are  significantly  smaller  than  their  counterparts  for  a  rigid  surface  or  rigid  bot¬ 
tom.  This  result  is  associated  with  the  fact  that  the  tweaking  term  in  the  group  velocity  is  positive  for 
free  or  negative  for  rigid  boundary  conditions. 

The  plots  of  group  velocity  vs  frequency  exhibited  the  familiar  behavior  of  going  to  zero  at  low 
frequencies  and  approaching  the  lower  of  the  boundary  sound  speeds  from  above  at  high  frequency. 

A  relative  maximum,  the  Airy  phase,  is  formed  for  intermediate  frequencies. 

The  profile  parameters  were  chosen  in  such  a  manner  that  the  physical  solutions  for  Case  A  and 
Case  B  should  be  almost  identical  when  the  surface  and  bottom  boundary  conditions  are  the  same. 
This  can  be  inferred  from  the  fact  that  here  Case  B  can  be  obtained  by  turning  the  Case  A  protile 
upside  down.  For  free-free  and  rigid-rigid  boundary  conditions,  the  Case  A  and  B  solutions  for  phase 
and  group  velocity  were  indeed  almost  identical,  even  though  the  intermediate  results  (such  as  eigen¬ 
values  or  canonical  group  velocity  factor)  were  markedly  different  for  Case  A  and  Case  B. 

The  eigenvalue  curve  for  mode  1  of  Case  A  and  Case  B  for  the  rigid-rigid  boundary  condition 
was  radically  different  from  that  of  any  other  mode  for  any  of  the  boundary  conditions.  However,  .*  • 
unusual  pattern  could  be  inferred  from  the  lattice  points  and  the  locus  of  the  higher-order  eigenval¬ 
ues.  Further  investigation  demonstrated  that  this  solution  was  not  only  mathematically  correct,  but  also 
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led  to  phase  and  group  velocities  which  made  sense  physically.  Here,  mode  1  corresponds  to  a  subset 
of  refracted  rays  which  do  not  reflect  from  both  surfaces.  In  contrast  to  all  other  modes,  there  is  no 
low-frequency  cutoff.  The  phase  and  group  velocities  go  to  a  finite,  nonzero  value  at  zero  frequency. 
The  square  of  this  value  is  the  harmonic  mean  square  of  the  sound  speeds  at  the  two  boundaries. 

This  value  is  slightly  less  than  the  arithmetic  mean  of  the  two  sound  speeds. 

The  bulk  of  the  eigenfunction  normalization  process  can  be  carried  out  in  terms  of  the  canonical 
eigenvalues  to  evaluate  the  canonical  eigenfunction  normalization  coefficient.  This  coefficient  applies 
to  the  entire  class  of  profiles  represented  by  the  eigenvalue  solution.  The  eigenfunction  normalization 
coefficient  for  a  specific  profile  of  the  class  is  obtained  by  a  simple  function  which  scales  the  canoni¬ 
cal  coefficient.  This  approach  is  more  straightforward  than  the  customary  approach. 

The  standard  approach  to  selecting  eigenfunction  coefficients,  before  normalization,  by  satisfying 
the  surface  boundary  condition  worked  well  for  the  Case  A  modes.  However,  the  Case  B  mode  nor¬ 
malization  became  very  large  for  higher  frequencies.  An  alternative  approach  was  implemented  in 
which  the  bottom  rather  than  surface  boundary  condition  was  satisfied.  Here,  the  Case  A  and  B  roles 
were  reversed,  with  a  satisfactory  result  for  the  Case  B  mode,  whereas  the  Case  A  mode  became  very 
large  for  higher  frequencies.  For  the  free-free  and  rigid-rigid  boundary  conditions,  the  eigenvalue 
normalization  coefficients  for  Case  A  (and  Case  B)  for  the  standard  approach  were  almost  identical  to 
those  for  Case  B  (and  Case  A),  even  though  intermediate  results  (such  as  the  canonical  normalization 
factor)  were  drastically  different. 

The  plots  of  the  eigenfunction  normalization  coefficients  vs  frequency  exhibit  a  characteristic  pat¬ 
tern.  The  values  are  always  positive.  The  coefficient  goes  to  zero  at  the  cutoff  frequencies  for  the 
free-rigid  and  free-free  boundary  conditions.  For  the  rigid-rigid  boundary  condition,  the  coefficient 
goes  to  for  the  cutoff  frequencies  for  modes  2  and  above.  Here,  Z2  is  the  duct  thickness.  At 

zero  frequency,  the  coefficient  goes  to  ji^Zi  for  mode  1  for  the  rigid-rigid  boundary  condition.  With 
the  exception  of  this  mode,  all  other  modes  for  all  boundary  conditions  form  a  relative  maximum  at 
intermediate  frequencies.  At  high  frequencies,  all  modes  asymptotically  go  to  zero.  For  the  rigid-rigid 
and  free-free  boundary  conditions  the  curves  for  the  various  modes  appear  to  merge  to  a  common 
curve  at  high  frequencies. 
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Appendix  A 

COMPARISON  OF  PHASE  VELOCITY  FOR  CASE  A  AND  CASE  B 

This  appendix  demonstrates  that  the  phase  velocity  vs  frequency  curves  are  almost  identical  for 
Case  A  and  Case  B  in  Fig.  9  and  10.  Let  the  variables  and  parameters  for  Case  A  be  given  by  x,  y, 
Cj,  C2,  and  y, .  Let  barred  values  represent  the  corresponding  values  for  Case  B.  The  Case  B  values 
are  related  to  Case  A  values  as  follows: 


*  =  y  ;  y  =  x  (a-i> 

-  Ci  ;  U2  =  Ci  ;  yi  =  -yt  (A-2) 

Equation  A-l  holds  because  the  curves  of  Fig.  3  and  4  are  symmetric  about  the  line  y  =  x.  Equation 
A-2  holds  because  of  our  choice  of  sample  profiles  for  Case  A  and  Case  B. 

We  now  compare  the  /  of  Eq.  14  of  Ref.  1  for  the  Case  A  values  with  the  /  of  Eq.  14.  Substi¬ 
tuting  Eq.  A-l  and  A-2  into  the  right  side  of  Eq.  14,  we  obtain  the  following  relationship: 

f  =  /(C,/C2)3  (A-3) 

If  we  now  substitute  Eq.  A-l  to  A-3  into  the  right  side  of  Eq.  17  of  Ref.  1,  we  obtain 

cp  =  c,[i  -  jrvY2/3\Yio\2/3yYy2  (A-4) 

A  comparison  of  Eq.  A-4  with  Eq.  30  yields  the  desired  result  that 

Cp  =  Cp  (A-5) 
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Appendix  B 

DUCT  CUTOFF  FREQUENCIES 

In  this  appendix,  we  obtain  expressions  for  the  duct  cutoff  frequencies  expressed  in  terms  of 
wavelength.  We  first  note  that 


/  =  N0C,Z-22  (B— 1) 

expresses  frequency,  where  N0  is  the  number  of  vertical  wavelengths  in  the  duct,  Ci  is  the  surface 
sound  speed,  and  Z2  is  the  duct  thickness.  Our  problem  is  to  obtain  expressions  for  N0  at  the  cutoff 
frequency  of  the  duct. 

From  Eq.  1  of  Ref.  1,  we  determine  that 

CiZ?  =  -2 y,(C?C;2  -  l)'1  (B-2) 

It  follows  then  from  Eq.  B-l  and  B-2  and  from  Eq.  14  of  Ref.  1  that 

No  =  ±  0  - Jf) 3/2  (C]C-22  -  1  )-1/2  (2*)-1  (B-3) 

The  plus  sign  holds  for  Case  B  and  the  minus  sign  for  Case  A.  Equation  B-3  holds  for  any  frequency. 

We  next  determine  the  conditions  corresponding  to  the  cutoff  frequencies.  The  cutoff  frequencies 
occur  when  Cp  of  Eq.  17  of  Ref.  1  goes  to  infinity.  This  occurs  when 

fn  =  *-2/3|y,|2'3x  (B-4) 

We  now  take  the  two-thirds  power  of  Eq.  14  of  Ref.  1,  equate  to  Eq.  B-4,  and  simplify  to  determine 
that 


jc  =  (y-x){C\C-2 2  -  I)'1  (B-5) 

We  take  the  one-half  power  of  Eq.  B-5  and  substitute  into  Eq.  B-3  to  obtain 

N0  =  ±  (y  -  x)x,/2(2n)~1  (B-6) 

where  the  plus  sign  holds  for  Case  E>  and  the  minus  sign  for  Case  A. 

Consider  now  the  lattice  points  of  Eq.  6  and  7  for  the  free-rigid  boundary  condition.  Let  the 
point  be  at  the  nt  zero  of  Eq.  6  and  the  n2  zero  of  Eq.  7.  We  consider  the  first  term  of  the  asymp¬ 
totic  expansion  of  Ai(-x)  as  given  by  section  10.4.60  of  Ref.  2.  If  we  set  this  term  to  zero,  we  find 
that 


or 


(2/3)x3/2  =  n(.nx  -  1/4) 


(B-7) 


x  =  (3  jt/2)2/3  n2/3(l  -  l/4ni)2/3 


(B-8) 
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We  next  consider  the  first  term  of  the  asymptotic  expansion  of  A/'(-y)  as  given  by  section  10.4.62  of 
Ref.  2.  The  counterparts  of  Eq.  B-7  and  B-8  are 

(2/3)y3/2  =  *(*2-3/4)  (B-9) 

and 

y  =  (3  */2)2/3  «2/J(l-3/4n2)2/3  (B-10) 

Consider  the  Case  A  modes.  Here 

n  =  *i  -  n2  +  1  (B-ll) 

where  n  is  the  mode  number.  Equation  B-ll  may  be  verified  by  noting  in  Fig.  2  the  relationship 
between  mode  number  and  lattice  point. 

We  next  use  the  binomial  expansion  on  the  last  factor  of  Eq.  B-8  to  obtain 

*  =  (3  */2)2/3  (ji\n  -  *i'/3/6)  (B-12) 

We  use  Eq.  B-ll  to  express  *2  >n  terms  of  *i  and  n  and  use  the  binomial  expansion  twice— once  on 
each  of  the  last  two  factors  of  Eq.  B-10  to  obtain 

y  =  (3  */2)2/3  [*?/3  +  2(l-n)*;,/3/3](l-nT‘/2)  (B-13) 

We  expand  Eq.  B-13  and  subtract  Eq.  B-12  to  obtain 

y-x  =  (3*/2)2/3*I1/3(l-2*)/3  (B-14) 

The  substitution  of  Eq.  B-12  and  B-14  into  Eq.  B-6  yields  the  final  result: 

N0  =  (2n  -  l)/4  (B-15) 

Consider  next  the  B  family  of  modes  for  the  free-rigid  condition.  The  demonstration  is  quite  simi¬ 
lar,  with  minor  changes,  as  will  now  be  noted.  The  counterpart  of  Eq.  B-ll  is 

n  =  *i  -  *2  (B-16) 

This  results  in  the  (1-n)  factor  in  Eq.  B-13  being  replaced  by  n.  This  in  turn  yields  the  negative  of 
Eq.  B-14,  and  the  result  of  Eq.  B-15  is  the  same  for  Case  B  as  for  Case  A. 

We  next  treat  the  free-free  boundary  condition  for  which  the  lattice  points  are  given  by  Eq.  13 
and  14  and  for  the  Case  A  modes.  Equation  B-12  applies  for  *  and  it  also  applies  for  y  with  n j 
replaced  by  n2.  The  counterpart  of  Eq.  B-ll  is 

*  =  *i  -  *2  (B-17) 


The  counterpart  of  Eq.  B-13  is 

y  =  (3  */2)2/3  [*?/3  -  2n«i1/3/3](l  ~n\'/6)  (B-18) 
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This  results  in  the  (1-2 n)  factor  of  Eq.  B-14  being  replaced  by  (-2 n).  The  counterpart  of  Eq.  B—  1 5 
then  becomes 


N0  =  n/2  (B— 19) 

We  next  treat  the  rigid-rigid  boundary  condition  for  which  the  lattice  points  are  given  by  Eq.  20 
and  21  and  for  the  Case  A  modes.  Equation  B-10  applies  for  y  and  it  also  applies  for  x,  with  n2 
replaced  by  nx.  The  counterpart  of  Eq.  B-12  is 

Jc  =  (3  n/2)111  (n?/J  -  n\u' /2)  (B-20) 

The  counterpart  of  Eq.  B—  1 1  is  the  same  as  Eq.  B-ll.  Thus,  Eq.  B-13  applies  for  this  case.  This 
results  in  the  (1  -  2 n)  factor  in  Eq.  B-14  being  replaced  by  2(1-  n).  The  counterpart  of  Eq.  B-15 
for  this  case  reduces  to 


N0  =  (n  -  l)/2  (B-21) 

For  the  Case  B  modes  for  the  free-free  and  rigid-rigid  boundaries,  the  signs  on  ni  and  n2  are 
the  reverse  of  that  in  Eq.  B-17  and  B-ll,  respectively.  This  results  in  a  sign  reversal  on  the  expres¬ 
sions  for  x  -  y.  Thus  Nc  for  the  Case  B  modes  reduces  to  the  same  expression  as  that  for  the  Case  A 
modes  for  both  sets  of  boundary  conditions. 

Before  closing  this  appendix,  a  few  comments  on  the  derivation  are  in  order.  The  reason  that  the 
analysis  is  so  intricate  and  involved  is  that  a  number  of  indeterminant  forms  must  be  dealt  with.  For 
example,  in  Eq.  14  of  Ref.  1,  we  must  deal  with  the  form  “>  -  »  as  both  x  and  y  increase  without 
bound.  Equation  B-14  and  its  counterparts  for  other  boundary  conditions  all  contain  the  factor  nj1/3- 
Thus,  we  see  that  as  ni  increases  without  bound,  the  frequency  goes  to  zero.  Since  the  lattice  points 
lie  on  the  curve,  it  follows  that  the  frequency  for  arbitrary  x  and  y  must  go  to  zero  for  x  and  y  with¬ 
out  bound.  We  chose  the  lattice  points  associated  with  the  zeros  of  Ai  and  A/'.  Alternatively,  we 
could  have  conducted  the  analysis  on  the  lattice  points  associated  with  the  zeros  of  Bi  and  Bi' .  This 
analysis  would  have  lead  to  the  same  result. 

Thus,  Eq.  B-14  demonstrates  our  statement  in  the  main  text  that  the  canonical  frequency  plots  go 
to  zero  asymptotically  as  x  becomes  large  without  bound. 

The  evaluation  of  Eq.  B-6  involves  another  indeterminant  form  of  the  type  zero  times  infinity.  In 
this  case,  y  -  x  depends  on  as  n\x,i,  whereas  the  dependence  of  x1/I  on  rt]  is  as  nj/3.  The  result 
is  that  Na  in  Eq.  B-15,  B-19,  or  B-21  is  a  finite  nonzero  constant. 

We  apparently  achieved  the  result  without  explicit]"  ’etting  n i  become  arbitrarily  large.  However, 
this  is  implicit  in  the  proof,  because  the  analysis  is  b  i  only  on  the  first  term  of  the  asymptotic 
expansions  of  Ai  and  Ai'  and  only  on  the  first  two  terms  of  the  various  binomial  expansions. 

In  closing,  Eq.  B-3  applies  for  any  x,y  solutions  of  the  eigenvalue  equation.  However,  Eq.  B-6 
does  not.  It  only  applies  at  the  cutoff  frequencies  as  given  by  the  condition  of  Eq.  B-5. 
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Appendix  C 

LIMITS  OF  THE  VARIOUS  NORMALIZATION  COMPONENTS 

Our  first  task  is  to  investigate  the  limits  of  the  various  expressions  for  Dc  given  in  the  main  text. 
Consider,  for  example,  G1  of  Eq.  44.  The  analysis  is  facilitated  by  a  simplification  of  Eq.  1.  We  solve 
the  eigenvalue  equation,  Eq.  2  and  5  for  Bi(-x),  substitute  for  Bi(-x)  in  Eq.  1  and  simplify.  The 
result  is 


(n  =  -  Ai(- x)/Ai'(- y)  (C-l) 

We  treat  the  G  functions  of  Eq.  46,  49,  52,  53,  and  54  in  the  same  manner,  using  the  appropri¬ 
ate  eigenvalue  equation  to  simplify  the  expression,  and  obtain,  respectively 


G2  =  n'  Ai(-x)/Ai'(-y)  (C-2) 

Gj  =  -  Ai’(- x)/Ai'(- y)  (C-3) 

Ga  -  -  ji~x  Ai’ (- y)/Ai{- x)  (C-4) 

G2  -  -jr”1  Ai’(-y)/Ai(-  x)  (C-5) 

and 

Gi  =  nx  Ai(-y)/Ai'(-  x)  (C-6) 


The  limits  of  Eq.  C-l  to  C-6  are  much  easier  to  evaluate  than  the  original  expressions  of  Eq.  1 
to  4.  Indeed,  in  the  computer  program,  Eq.  C-l  to  C-5  were  used  to  reduce  the  effect  of  numerical 
roundoff. 

Consider  now  the  limits  for  large  positive  x  and  y.  Here  y  -*  x  asymptotically  and  the  Airy  func¬ 
tion  ratios  in  Eq.  C-2,  C-3,  C-5,  and  C-6  go  to  1.  From  Sections  10.4.59  and  10.4.61  of  Ref.  2, 
the  Airy  function  ratio  in  Eq.  C-l  goes  to  -  (xy)"1/4  and  to  -  (xy)1/4  in  Eq.  C-4.  It  is  evident  then 
that  the  second  terms  of  Eq.  44  and  46  go  to  n~2,  the  second  term  of  Eq.  49  goes  to  y: r'2,  the  sec¬ 
ond  term  of  Eq.  52  and  54  goes  to  -xn'2,  and  the  second  term  of  Eq.  53  goes  to  -n'2.  Thus,  Eq. 
44,  46,  and  53  reduce  to 


n 


-2 


ji-2  =  0 


(C-7) 


whereas  Eq.  49,  52,  and  54  reduce  to 


x{n2  -  n~2)  =  0  (C-8) 

It  is  of  interest  that  as  x  becomes  large  positive,  Fig.  17,  18,  and  24  approach  zero  much  faster 
than  do  Fig.  19,  23,  and  25.  The  reason  is  that  the  former  three  correspond  to  Eq.  C-7.  The  latter 
three  correspond  to  Eq.  C-8,  in  which  the  difference  of  functions  not  only  must  cancel  each  other  as 
in  Eq.  C-7,  but  also  must  drag  down  the  multiplier  x  of  Eq.  C-8. 

We  next  examine  the  limits  for  large  negative  x  or  y.  Recall  that  for  the  Case  A  modes,  x  goes  to 
the  zeros  of  Ai(-x)  or  Ai'  (-x)  as  y  goes  to  minus  infinity;  whereas  for  the  Case  B  modes,  y  goes  to 
the  zeros  of  Ai(-y)  or  Ai'(-y)  as  x  goes  to  minus  infinity.  Thus,  Eq.  C-l  to  C-6  are  all  indeterminate 
forms  of  the  type  0/0. 
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We  have  attempted  to  evaluate  these  forms  by  applying  L’  Hospital’s  rule  without  success.  This 
procedure  is  complicated,  because  it  involves  dy/dx,  and  one  differentiation  did  not  provide  the 
answer.  However,  we  will  make  one  relatively  simple  assumption  and  then  demonstrate  that  all  of  the 
numerical  results  are  compatible  with  this  assumption.  The  assumption  is  that  a  finite  zero  of  an  Airy 
function  goes  to  zero  faster  than  does  the  zero  at  x  or  y  equal  to  minus  infinity. 

This  assumption  leads  to  zero  and  infinity  as  the  limits  of  Eq.  C-l  to  C-3  for  Case  A  and  Case  B 
modes,  respectively.  This  is  compatible  with  the  results  of  Fig.  17  to  19,  in  which  the  Case  A  modes 
go  to  finite  limits  and  the  Case  B  modes  appear  to  diverge.  Eq.  C-4  to  C-6  involve  the  reciprocals  of 
the  right  sides  of  Eq.  C-l  to  C-3,  respectively.  Thus,  the  roles  are  reversed,  and  the  assumption  leads 
in  Fig.  23  to  25  to  finite  limits  for  Case  B  and  apparently  divergent  values  for  Case  A. 

We  next  examine  the  limits  of  Dn  for  high  frequency.  This  involves  the  limits  of  Dc  as  x  or  y 
goes  to  minus  infinity,  as  just  discussed.  For  those  cases  in  .vhich  Dc  has  a  finite  limit,  it  is  clear  from 
Eq.  89  and  from  Eq.  4  of  Ref.  1  that  D„  goes  to  zero  as  suggested  by  the  A  modes  in  Fig.  20  to  22 
and  the  B  modes  in  Fig.  26  to  28.  Apparently  the  j~2/y  dependence  in  a^1  is  not  strong  enough  to 
drag  down  a  divergent  Dc .  Thus,  D„  apparently  diverges  for  the  B  modes  in  Fig.  20  to  22  and  for 
the  A  modes  in  Fig.  26  to  28. 

Our  final  task  is  to  examine  the  limits  of  Dn  as  the  frequency  decreases  to  low  values.  From  Eq. 
14  of  Ref  1  and  Eq.  B-2,  it  follows  that 

//3  =  -  (y-jc)C,  Z;1  yi1/3  (C— 9) 

Substitution  into  Eq.  4  of  Ref.  1  leads  to 


ai  =  (y-x)z~2x  (C- 10) 

Thus 

Dn  =  Z2  Dc/(y-x)  (C-l  1) 

Dc  is  of  the  form  n-2  (zero)  for  the  free-rigid  and  free-free  boundary  conditions  and  ji~2(y-x) 
for  the  rigid-rigid  boundary  conditions. 

Consider  now  the  case  of  the  modes  1A  and  IB  for  the  rigid-rigid  boundary  condition.  Here, 
x  =  -y  at  low  frequencies  and  x  approaches  zero.  The  result  then  appears  to  reduce  to 

Dn  =  Z2  a1  (C-l 2) 


which  is  the  limit  shown  in  Fig.  22. 

However,  for  higher-order  modes,  x  -*  y  asymptotically.  A  similar  analysis  appears  to  give  the 
same  answer  as  Eq.  C-12.  However,  the  numerical  lir  it  shown  in  Fig.  22  is  one  half  of  Eq.  C-12. 
We  puzzled  over  this  result  until  we  recognized  that  at  duct  cutoff,  Eq.  B-4  holds.  If  we  use  Eq.  B-4 
to  evaluate  a  i ,  we  obtain 


a\  =  -2  y,  jr/Ci 


(C-l  3) 


For  the  free-rigid  and  free-free  boundary  conditions 

D„  =  -Ci  zr'!  {zero)/ 2  yx  x  (C— 14) 
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As  x  goes  to  infinity,  Eq.  C-14  goes  to  zero,  in  agreement  with  Fig.  20  and  21.  For  the  rigid-rigid 
boundary  condition 


D„  =  -Ci  y,  x  (C-15) 

Here,  as  y  approaches  x,  the  limit  of  Eq.  C-15  appears  to  go  to  zero,  which  does  not  agree  with  Fig. 
22. 


The  problem  here  is  that  Eq.  C-15  does  not  take  into  account  that  Eq.  B-6,  as  well  as  Eq.  B 
must  hold.  If  for  the  rigid-rigid  condition,  the  limit  of  Dc  as  y  approaches  x  asymptotically  is 

Dc  =  sr2(y-x)/ 2  (C-16) 

and  the  result  of  Fig.  22  is  obtained.  The  proof  of  Eq.  C-16  appears  to  require  an  analysis  of  lattice 
points  similar  to  that  of  Appendix  B  and  appears  too  difficult  to  attempt  here. 

In  the  discussion  following  Eq.  B-4,  it  was  noted  that  the  standing  wave  for  the  rigid-rigid  bottom 
at  cutoff  was  approximated  by  a  cosine  function.  If  one  normalizes  this  cosine  function,  one  obtains, 
Z2/2  as  the  normalization.  The  Airy  functions  involve  a  factor  of  jt~l  times  the  trigonometric  func¬ 
tions.  This  leads  to  the  desired  result: 


Dn  =  Z2  n1 /I  (C-17) 

Although  the  various  limits  of  Dc  and  D„  are  of  interest  as  a  mathematical  challenge,  our  pri¬ 
mary  reason  for  the  analysis  in  this  appendix  was  a  verification  of  the  numerical  results.  We  wish  to 
guard  against  possible  error  in  the  numerical  implementation  of  the  theory.  Whereas  our  mathematical 
analysis  leaves  some  significant  unresolved  gaps,  there  are  enough  points  of  agreement  to  make  us  rea¬ 
sonably  certain  that  the  numerical  results  are  not  in  error. 
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